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Abstract 

Classically it was held that solutions to deterministic partial differential equations 
(i.e. ones with smooth coefficents and boundary data) could become random only 
through one mechanism, namely by the activation of more and more of the infinite 
number of degrees of freedom that are available to such a system. It is only recently 
that researchers have come to suspect that many infinite dimensional nonlinear systems 
may in fact possess finite dimensional chaotic attractors. Lyapunov exponents provide 
a tool for probing the nature of these attractors. In this paper we examine how these 
exponents might be measured for infinite dimensional systems. 
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1 Introduction 


An attracting set for a dynamical system is a region in phase space which ‘attracts’ nearby 
initial conditions. Any orbit started in the neighbourhood of such a set will evolve towards 
it and not leave it thereafter. Attracting sets for dissipative systems have dimensions which 
are less than those of the phase space as a whole and as they eventually trap all initial 
conditions it is their character which governs the long term, asymptotic behaviour of the 
system. It is clearly in our interest to characterize these sets as closely as we can. Attracting 
fixed points and periodic orbits provide two examples of ‘well behaved’ (and consequently 
dull) attracting sets. Recently however, examples have been found of some remarkable 
attractors. These strange attractors are characterized by the fact that orbits in them, 
which at some time lie infinitely close together, diverge from each other at an exponential 
rate and become uncorrelated in a finite time. Such behaviour, termed sensitive dependence 
on initial conditions, must be present if the system is to be considered chaotic. 

Many finite dimensional flows are now known to have this sensitive dependence on 
initial conditions. Perhaps the best known of these is the flow of the Lorenz [1] set of 
three nonlinear ordinary differential equations (for appropriate parameter values) . All the 
initial conditions quickly settle onto the attracting set but if one evolves an infinitesimal line 
segment of initial conditions these quickly smear out over the entire attractor. The initial 
conditions for any experiment, whether it be on the computer or in the laboratory, can 
only be specified to a finite precision and if the attractor for the system ha s the exponential 
orbital divergence property, states which are initially so close that we cannot resolve their 
differences will, in a finite time, give rise to quite different behaviours. 

It is important to note that while we speak of exponential divergences the motions may 
in fact be taking place in a bounded domain. Thus two nearby orbits will diverge but may 
at a later time lie once more quite close together and never get infinitely far apart. The 
important notion here is one of the loss of predictability and it is this notion which we seek 
to quantify by means of the Lyapunov spectrum of the system. The rest of this paper is 
given over to a discussion of the Lyapunov exponents, their relevance to this problem, and 
in particular, their relevance to chaos in infinite dimensional dynamical systems. 
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2 Lyapunov Exponents 


Lyapunov exponents (also called characteristic exponents) essentially measure the mean 
exponential rate of divergence of nearby orbits of a dynamical system. Let us illustrate the 
basic idea by considering a discrete time dynamical system, the map 

x n + 1 = /(^n)* (1) 

Orbits are generated by iterating on an initial seed xo- Consider two nearby initial seeds 
xo and Xq separated by an infinitesimal distance 8. Assuming exponential divergence at a 
rate A (which may in general depend on xo ) we have after N iterations, 

5e A(.o)I» = !,"(*„) _ /»(4)|, ( 2 ) 

whence 

M*o) « lo 8 Df N (x 0 ), (3) 

where D represents the derivative operator. Applying the chain rule we get that 

j N—i 

A M « T7 ^2 log Df(x n ). (4) 

n= 0 

This suggests that we define the Lyapunov exponent for the map as follows, 

A(*o) = ^ log Df(x n ). (5) 

In this fashion we examine the linearization of the map along the trajectory generated 
by the seed xo- The Lyapunov exponent A(xo) measures the mean exponential expan- 
sion/contraction rate along this orbit. Oseldecs [2] theorem can be invoked to address the 
question of the existence of the limit in equation (5). The interested reader is referred 
to [3] and the references contained therein for further details. 

Next we consider a system of n ordinary differential equations, 

i = /(*), (6) 
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where x = {xi , . . . , x n } . 

Assuming the system is dissipative any initial n dimensional volume in phase space 
will be contracted by the flow of (6). We will consider the evolution of spheres of initial 
conditions centered about some orbit which we will term the fiducial (faithful) trajectory. 
Exponents are measured then for a particular fiducial orbit. However there is usually an 
assumption made about the ergodicity of the system. Roughly speaking we assume that we 
would get the same values for the exponents no matter which initial condition was chosen 
for the fiducial trajectory. For this to be true the system must have an indecomposable 
attractor which is visited almost everywhere by almost all orbits. That is, it must not 
be possible to decompose the attractor into distinct parts which might well have different 
properties. Moreover, almost all orbits should fill out the attractor so that the calculations 
(at least if performed over a long enough time span) will not depend on the particular piece 
of the attractor the fiducial trajectory visits. 

A sphere of initial conditions will be deformed into an ellipsoid whose volume (which 
is proportional to the product of its semiaxes) will be smaller than the original. This does 
not preclude the possibility that some of the semiaxes are in fact larger than the radius of 
the original sphere. Initial conditions lying along such an axis will diverge from each other. 
Of course as the system evolves in time the semiaxes of the ellipsoid will in general not 
point in some constant direction but they will rotate as well as translate. This means that 
the detailed structure of the orbital divergences will be very complicated but the Lyapunov 
exponents shall be defined to give us a method of quantifying these divergences without 
considering all the complications 

The largest Lyapunov exponent will measure the mean exponential rate of maximal 
orbital divergence along a trajectory. To put this another way we think of ourselves moving 
along our fiducial trajectory making measurements at each point. We have at our disposal 
an elastic (stretchable) meter stick, one end of which is tied to the trajectory with the other 
end being left free to evolve under the influence of the differential equation. The free end 
of the stick will tend to line up with the direction of maximal expansion of the system and 
it is the resultant change in its size that we measure. The first Lyapunov exponent is then 
defined as the mean value of the logarithm of our measurements along the orbit. 

The higher order exponents are defined by considering the evolution not of lines but 
rather of volumes of initial conditions. The first two exponents will then measure the mean 
value of the maximal stretching of a disc of initial conditions. The circular disc deforms 
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into an elliptical one whose area is proportional to the product of its two semiaxes. The 
first two exponents refer to the mean exponential rate of change of the lengths of these axes 
along the fiducial trajectory. 

Consider a general m-dimesional ellipsoid with time dependent semiaxes ai(t), . . . , a m{t)- 
The volume of the ellipsoid is proportional to w m (t) where w m (t) = aci(t)a2[t) . . . a m (t). The 
Lyapunov exponent associated with a,-(t) is then defined to be 


A, 


lim 

f— *oo 


log <*«(*) 


(7) 


provided the limit exists. 

The rate of growth of a surface element is given by the sum of the largest two exponents, 

Ai + A2, while the rate of growth of a fc-volume element is given by the sum, Ai + 1- A*. 

We note that for a dissipative system the action of the flow is to contract the phase space 
as a whole and therefore the sum of all the exponents will be negative. 


In the infinite dimensional case we have an equation of the form, 

u = N(u), (8) 

where N is some nonlinear partial differential operator and u is a function of many variables. 
We will assume that u belongs to some normed space with a basis and for simplicity let us 
assume that this basis defines the boundary data for u (i.e. each basis function individually 
satisfies the imposed bundary conditions). We now consider the problem we had before, 
namely, whether nearby initial conditions evolving under the influence of (8) diverge or 
contract. 

Distances are measured in the appropriate norm and we emphasise that any contraction 
or expansion is taking place in the function space and is generally not easy to interpret 
in physical space. Expansion can be thought of as a modal amplification with the la- 
belling of the amplified modes defining the direction of expansion in the function space. 
Whereas in the finite dimensional system we might have spoken of solutions diverging in 
the ‘23’ direction we might now speak for example of divergences in the ‘sin 2x3’ direction. 

Consider two nearby initial distributions, uo(z) and Uq(x) for the unknown in equation 
(8) evolving into the solutions u(x, t) and u f (x,t) respectively. Define S(x,t) to be the 
difference, 

S(x y t) = u(z, t) - u'(x,t). (9) 
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Divergence of the solutions in some direction will manifest itself as an amplification of the 
corresponding modes in the basis function expansion of 8(x,t). As in the finite dimensional 
case the direction of divergence may change in a complicated manner with time. That is, 
different modes may get amplified at different times. The Lyapunov spectrum will again be 
defined so as to avoid this complication of orientation and it will merely track the norm of 
the divergences. 

In fact the exponents are defined by considering the evolution of volumes of initial 
conditions just as was done in the finite dimensional case. A k-sphere of initial states in 
function space is realized by taking any k basis functions as k separate initial conditions 
(assuming the basis functions are orthonormal). We can watch such a sphere evolve and 
deform under the influence of (8) . Even though we initially take rather simple, single mode 
states the nonlinearities in the governing equation will generate new modes in the solutions 
as time evolves. In our picture the appearance of these new modes correspond to rotations 
of the initial sphere while changes in amplitudes of the modes correspond to deformations 
of the sphere into an ellipsoid. 

The volume of the ellipsoid is again proportional to the product of its semiaxes with 
the semiaxes now being thought of as modal amplitudes. The Lyapunov exponents are 
defined in precisely the same manner as was done above in (7). We will discuss how such 
a calculation might be carried out in practise in section 5. 

3 Dimension 

It was long thought that solutions to deterministic partial differential equations could be- 
come random only through one mechanism, namely by the activation of more and more 
of the infinite number of degrees of freedom that are available to such a system. Recently 
however, researchers have come to suspect that many infinite dimensional systems may in 
fact possess finite dimensional chaotic attractors. For a finite set of ordinary differential 
equations, theoretically at least, one can measure all the characteristic exponents for the 
system and if their sum is negative then regions in phase space are contracting. Then not 
all the possible degrees of freedom are active and we might say that the attractor is of a 
lower dimension than the phase space as a whole. 

The descent from an infinite dimensional phase space to a finite dimensional attractor 
is considerably more difficult to make. Early work in this area was done by Foias, Manley, 
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Temam and Treve [4]. The important question is when will a finite approximation to the full 
solution of a partial differential equation give qualitatively correct results. The approxima- 
tions that one has in mind here consist of finite sums of appropriate basis functions. Such 
a numerical methodology is of course well known and goes under the general appellation of 
spectral methods (see for example [6]). By qualitative agreement we mean at what level of 
truncation will the finite approximation and the full solution display the same behaviour 
with regard to such properties as stability, periodicity, quasi-periodicity etc. 

The approach taken is to try and show the existence of a finite dimensional Lipschitz 
manifold, an inertial manifold , for the equation in question. The distinguishing properties 
of this manifold are: 

1. It is invariant under the flow of the system. 

2. All the solutions to the equation converge exponentially to the manifold no matter 
where they are initiated in phase space. 

Systems which have been shown to possess an inertial manifold include the one dimensional 
Kuramoto-Sivashinsky equation which models turbulent fluid interfaces (for an extensive 
treatment of this equation see [7]). This partial differential equation is then completely 
equivalent to a finite system of ordinary differential equations. 

We note in passing that the proof of existence of the inertial manifold is reminiscent of the 
proof of the existence of a center manifold in the finite dimensional context. However, where 
the center manifold is a local construction, tied to a particular fixed point of the ordinary 
differential equation, the inertial manifold is a global object for the partial differential 
equation. It is the local nature of the center manifold which make it useful in calculations, 
as, close to the fixed point, it can be approximated by polynomial curves (a Taylor series 
approach). With such a specific representation of the lower dimensional surface at hand it 
is a relatively simple matter to unfold the system in this neighbourhood of interest and to 
determine the nature of the flow there (see for example [5]). The inertial manifold on the 
other hand is a global object which does not provide any clues as to which are the points 
of interest that one should try to expand about. As yet there is no systematic procedure 
for making use of the inertial manifold beyond saying that if it exists the system is in some 
strict sense equivalent to a finite dimensional system. 

Any attracting sets will be subsets of the inertial manifold and it is these sets which 
are all important in the asymptotic limit. Unlike the manifold these sets may have no 
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smoothness properties to speak of. There are generalizations of the concept of dimension 
that may be usefully invoked when we discuss the extent and properties of these attractors. 
These generalized dimensions can be thought of as quantifying the amount of information 
that must be provided to specify a point on the attractor to within given error bounds. 
Extensive discussions of dimension can be found in Mandelbrot [8]. 

The many definitions of dimension fall for the most part into two camps, the probabilistic 
and non- probabilistic definitions (for a longer discussion the reader is referred to [9]). The 
former take into account the frequency with which the various parts of the set are visited 
by a typical trajectory whereas the latter merely measure metric properties of the set and 
ignore the dynamical features. A simple and appealing notion of a dimension which is 
non-dynamical in nature is that of the capacity of a set X which is defined as 


dc = 


log^(e) 

~o log(l/e) * 


( 10 ) 


where N(e) is the number of balls of radii < e needed to cover X. Intuitively we are 
defining dc so that the volume of the set is proportional to Ne d °. For simple sets such 
as points, lines and areas it is easy to verify that dc yields the usual results of 0, 1 and 2 
respectively. However for ‘strange’ sets dc will yield non-integer dimensions. For instance 
for the well known case of the middle thirds Cantor set in one dimension, dc = 0.630 . . ., a 
result that fits in with the notion that this set is something less than a line and something 
more than a point. 

To include the dynamics and the frequency of visitation in our discussion of dimension it 
is useful to define a quantity called the natural measure. Consider a ball B in the attractor 
and an initial condition x in the basin of attraction. Define p(x, B) as the fraction of time 
that a trajectory started from x spends in B. If p(x, B) is the same number for almost all 
x in the basin of attraction then we denote it by p{B) and call p the natural measure of the 
attractor. 


By using such a measure to weight the balls in (10) appropriately we can ask what the 
capacity of the most probable part of the attractor is. In other words we can require that the 
N(e) balls used in (10) be chosen so that their natural measure adds up to some fraction, 6 
of the measure for the full attractor (which by definition is 1). The capacity measured this 
way is called the 0-capacity and is denoted dc{9). Clearly dc(l) = dc, however the value 
of dc{6) for 0 < 1 may be different from its value at 9 = 1 . Experimentally it has been 
found that for many systems dc(0) is independent of 0 and that dc(0) < dc for 6 < 1. 
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4 Using Exponents to Estimate Dimension 


Recently there has been work done which seeks to delineate the connection between Lya- 
punov exponents and the dimension of attractors. A heuristic argument based on the case 
of an attractor having just one positive and one negative exponent, and A 2 respectively, 
may shed some light on a possible connection. Begin by taking a covering of the attractor 
with A(c) discs, each of radius c. For small enough € the action of the flow is roughly linear 
over the interior of the discs for some finite time r. In this time the discs will have deformed 
into ellipses with areas proportional to A finer covering of the attractor with 

discs of the smaller radius ee* 3T can also be considered and it will take roughly 
such discs to cover one of these ellipses. If we use the assumption that each evolved ellipse 
has the same area, we arrive at the estimate 


N(ee X2T ) = e (Al_A2 ) r iV(e). (11) 

The definition of the capacity suggests that N (e) « fc(l/e) d ° and if we substitute this 
relation into the last equation we can solve for dc 


‘ fe = 1_ XT 

Thus motivated, we define the Lyapunov dimension, di for this case by 


di, — 1 


Ai 

A2 


The generalization to the higher dimensions is given by 


di — m- 


E£iA, ; 

A«i+i 


( 12 ) 


(13) 


(14) 


where m is the largest integer for which A,- is positive. 

This dimensional quantity is due to Kaplan and Yorke [10]. While we have drawn a 
connection between the Lyapunov dimension and the capacity, it must be noted that the 
Lyapunov exponents are mean quantities and, as such, are affected by the natural measure 
of the attractor. It has been conjectured [9], that for a typical attractor di = dc{9) with 
9 < 1. 


A useful result linking the capacity and the Lyapunov exponents has been proved by 
Constantin and Foias [11], This result is given in terms of the uniform Lyapunov exponents, 
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/i» which are defined iteratively by the equation 


Hi = lim (log -^|r ] ft. 

The relation puts an upper bound on the capacity as follows 

d c < max (/ + fi =1 ) , 

|A m+ i| J 

where m is such that Hi < 0 and strictly speaking we must take the supremum of 

each quantity over all the trajectories in the attractor for equation (16) to hold. 

This result can be used in conjunction with the knowledge of a large enough portion of 
the Lyapunov spectrum to find bounds on the dimension of the attractor. The best bound 
will be obtained by taking m to be the last integer for which YliLi Hi > 0. Thus we should 
like to be able to calculate Lyapunov exponents for higher and higher dimensional volumes 
until we find a volume that is contracted by the flow of the system and then we can use (16) 
to put a bound on the capacity of the attractor. 

5 Measuring the Exponents 

It is possible in theory to measure the first m Lyapunov exponents by measuring the rate of 
separation of m nearby trajectories. However there are many numerical problems associated 
with this method. The exponents supposedly probe the local structure of the attractor, 
being defined in terms of the mean values of some small growths in distances. Many chaotic 
attractors have embedded in them strange folds and twists, for example the ‘double ear’ 
structure of the Lorenz attractor. In the Lorenz case, while monitoring nearby orbits, we 
may at some time see a separation where one of the trajectories begins to traverse the right 
ear of the attractor while the others continue to travel around the left ear. At the time 
of separation we will witness a sudden jump in inter- trajectory distances. This jump is 
indicative of the global geometric properties of the attractor and it is not what we wish to 
quantify with the Lyapunov exponents. 

A method which was discovered independently by Bennetin et al. [12] and Shimada 
and Nagashima [13] gets around such difficulties. This method advocates following just 
one trajectory in phase space (termed the fiducial trajectory earlier) and simultaneously 


( 15 ) 


(16) 
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following the trajectories of m independent vectors in the tangent space to the fiducial 
trajectory. These evolve under the influence of the associated linearized equations. In 
this fashion we are looking at the evolution of an infinitesimal sphere of initial condition 
centered at all times on the fiducial trajectory. If the fiducial trajectory takes a particular 
path around some part of the attractor it will drag the other m linearized trajectories along 
that path with it. In this way we can be sure that we continue to measure local properties 
of the attractor at all times. 

If the equation for the fiducial trajectory, u(t) is given by (8) then the equation for a 
slave, tangent space, orbit, v(t) is 

v = DN(u)v, (17) 

where DN(u) represents the linearization of the nonlinear operator, N, about u(t) To 
calculate m exponents we solve (8) for u(t) and simultaneously solve (17) for m different 
initial conditions, vi 0 , • • • > t>n»o- Tracking the divergences in vi(t), . . . , v m (f) gives us the 
exponents. 

However, there are still numerical difficulties. The exponential growths in certain com- 
ponents of each of the tangent vectors may result in a decrease in the angle subtended by 
the vectors. Rather quickly on any finite precision computer these vectors will appear to 
become aligned and consequently are no longer independent. When this happens all of the 
vectors monitored will give the largest Lyapunov exponent as the rate of its growth. The 
familiar Gram-Schmidt orthogonalization procedure can be profitably employed to avoid 
this duplicity. 

The modified approach is to once again solve (8) for u(t) and to simultaneously solve 
m copies of the linearized equations, (17) for vj(t), . . . , v m (t). The initial conditions, 
vi 0 , . . . , Vmo ar e taken to be orthonormal. Not only that but the solutions vi (t) , . . . , v m (t) 
are not let to evolve as they wish but are reorthonormalized every so often by means of the 
Gram-Schmidt procedure. 

Under the flow the first vector, «i(t) tends to line up locally with the direction cur- 
rently associated with the largest Lyapunov exponent. The Gram-Schmidt procedure merely 
shrinks this vector and leaves its orientation unaltered. The surface element defined by the 
first pair of vectors, vj(t) and V 2 (t) aligns itself along the plane defined by the eigen-directions 
associated with the largest two exponents, Aj and A 2 . The Gram-Schmidt procedure alters 
both the magnitude and the orientation of V 2 (i). However, the effect of this is to merely 
alter the area of the surface element leaving it aligned as it was before. Similarily, the 


10 


higher dimensional volumes will not have their orientation altered by the procedure though 
the defining vectors will be rotated. We remind the reader that it is the volume elements 
that are important for the exponent calculations. 

The effect of the Gram-Schmidt procedure is twofold. It prevents the tangent vectors, 
which are supposedly defining volumes, from degenerating and collapsing onto each other 
and it also renormalizes exponentially growing numbers (vector norms), thus inhibiting 
overflow problems on the computer. Moreover, the norms calculated as part of the procedure 
are precisely the quantities we wish to monitor in our exponent calculations. Wolf et al. [14] 
can be consulted for the details of how this method can be used for finite dimensional 
systems. 

Here we are interested in infinite dimensional systems and as an example we will now 
outline how the method can be used to find the exponents associated with the Kuramoto- 
Sivashinsky equation: 

u t = -Au xxxx - au xx - ^(u*) 1 2 . (18) 

where a is a real parameter. We assume that u satisfies perodic boundary conditions in x 

u(x + 2?r, t) = u(x, 0). (19) 

A suitable, high resolution, numerical technique for solving (18) is the Fourier collocation 
method. We will not go into the details of such a scheme here but merely assume that some 
such algorithm has been implemented on the computer and yields acceptably accurate, 
approximations, u(x,t), to the true trajectories. 

The associated linearized equation for a slave trajectory t>(x, t) is 

Vf = 4Vxzxz GWxx OLMx^x* (20) 

and again v(x, t) is 2n periodic in x. We assume that a suitable means of solving (20) is 
available and that it too yields approximate solutions v(x,t). The Gram-Schmidt procedure 
must also be implemented numerically; a simple matter if the Fourier collocation method 
is being used to solve for the trajectories. 

With these numerical tools in place we proceed to calculate m Lyapunov exponents as 
follows. 

1. Starting from an arbitrary initial condition we evolve u(x, t) according to equation (18) 

for some time, T. The value, T, must be picked (by trial and error) to be large enough 
that we can feel confident that u(x,T) is close to the attractor. 
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2. We then start to solve (20) for the slave trajectories. Take orthonormal initial con- 
ditions for these, for example, vjfe 0 (x) = sin(fca;), k = 1,.. . ,m. Begin by setting the 
t7jk(a:,0) equal to the chosen initial conditions. 

3. Update the solution u(x,t) by one further time step using equation (18). 

4. Update each slave trajectory, v k (x, t), k = 1, . . . , m for one time step using the values 
of u(x,t) calculated previously, in place of u(z,t) in the governing equation (20). 

5. Perform the Gram-Schmidt procedure on the v’s. Keep cumulative records of the 
logarithms of each of the new norms. The current estimates of the Lyapunov exponents 
are given by these records divided by the time elapsed since we started to solve 
equation (20). 

6. Go to step 3. 

Of course there is no way of telling when the fiducial trajectory will have settled onto the 
attacting set. Nor indeed is there any guarantee that there is just a single attractor. What 
we can do, is to experiment with various initial conditions for the fiducial trajectory and 
also try various values for the settle time, T. If the calculated exponents prove robust to 
such experimentation we can hope that they reflect the structure of the existing attractor. 
In any case it is our intention here to just present a tool for such investigations. We leave 
it for another time to consider the results that can and cannot be extracted with this tool. 

We would like to conclude by considering the consequences of the existence of a positive 
Lyapunov exponent for an infinite dimensional system. There is a famous criterion of 
Hadamard for the well posedness of a partial differential equation. This essentially states 
that in order to be physically meaningful, equations such as (8) should be stable in the sense 
that small changes in the initial data should produce only small changes in the solutions. 
However it is now felt that there are many physical systems which do not in themselves 
obey this criterion. If the model equation for the system has sensitive dependence on initial 
conditions this may well be indicative of a physically significant phenonema. 

Nor does the the existence of some postive Lyapunov exponents doom the numericist 
to failure. While we are constrained at all times to a finite precision calculation and con- 
sequently cannot be certain of precisely repeating a given experiment (say of repeating a 
given computer run on a different machine ), we can hope to calculate accurately, mean 
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quantities for the system either by averaging over many runs or over large times for one 
run. Indeed in any practical sense such quant it es are bound to be of more value than a 
knowledge of the detailed structure of a chaotic attractor which will be extremely complex 
and difficult to interpret. 
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